#packages
library(haven)
library(lattice)
library(misty)
library(Matrix)
library(lme4)
library(msm)
library(sandwich)
library(sjPlot)
library(lmerTest)
library(backports)     
library(carData)
library(effects)       
library(ggplot2)       
library(interactions)  
library(psych)         
library(plyr)
library(ggpubr)

#################################################################################################################
#Prepare data for base model
########################################################################################################################

dataset_rmssdbase <- dataset[ , c("id", "stressmarker", "resting_rmssd", "reactivity_rmssd", "resting_met", "reactivity_met", "conttime")]
dataset_rmssdbase <- dataset_rmssdbase[complete.cases(dataset_rmssdbase),]

#RMSSD reactivity
dataset_rmssdbase$rmssd_reactivity <- (dataset_rmssdbase$reactivity_rmssd - dataset_rmssdbase$resting_rmssd)

#Log transform RMSSD parameters
#RMSSD reactivity
dataset_rmssdbase$reactivity_rmssd_log <- log(dataset_rmssdbase$reactivity_rmssd)
#resting RMSSD 
dataset_rmssdbase$resting_rmssd_log <- log(dataset_rmssdbase$resting_rmssd)

#Center independent variables
#RMSSD resting
dataset_rmssdbase$resting_rmssd_log_cent <- center(dataset_rmssdbase$resting_rmssd_log, type = c("CWC"), cluster = dataset_rmssdbase$id, value = NULL, as.na = NULL,
       check = TRUE)
#Reactivity MET
dataset_rmssdbase$reactivity_met_cent <- scale(dataset_rmssdbase$reactivity_met, scale = FALSE)
#Resting MET
dataset_rmssdbase$resting_met_cent <- scale(dataset_rmssdbase$resting_met, scale = FALSE)
#Conttime_hour
dataset_rmssdbase$conttime_cent <- scale(dataset_rmssdbase$conttime, scale = FALSE)

#rescale
dataset_rmssdbase$conttime_cent <- dataset_rmssdbase$conttime_cent/24

#outliers
#Reactivity RMSSD
dataset_rmssdbase$rmssd_reactivity[which(dataset_rmssdbase$rmssd_reactivity < (mean(dataset_rmssdbase$rmssd_reactivity) - 3*sd(dataset_rmssdbase$rmssd_reactivity)))] <- (mean(dataset_rmssdbase$rmssd_reactivity) - 3*sd(dataset_rmssdbase$rmssd_reactivity))
dataset_rmssdbase$rmssd_reactivity[which(dataset_rmssdbase$rmssd_reactivity > (mean(dataset_rmssdbase$rmssd_reactivity) + 3*sd(dataset_rmssdbase$rmssd_reactivity)))] <- (mean(dataset_rmssdbase$rmssd_reactivity) + 3*sd(dataset_rmssdbase$rmssd_reactivity))
#RMSSD resting
dataset_rmssdbase$resting_rmssd_log_cent[which(dataset_rmssdbase$resting_rmssd_log_cent < (mean(dataset_rmssdbase$resting_rmssd_log_cent) - 3*sd(dataset_rmssdbase$resting_rmssd_log_cent)))] <- (mean(dataset_rmssdbase$resting_rmssd_log_cent) - 3*sd(dataset_rmssdbase$resting_rmssd_log_cent))
dataset_rmssdbase$resting_rmssd_log_cent[which(dataset_rmssdbase$resting_rmssd_log_cent > (mean(dataset_rmssdbase$resting_rmssd_log_cent) + 3*sd(dataset_rmssdbase$resting_rmssd_log_cent)))] <- (mean(dataset_rmssdbase$resting_rmssd_log_cent) + 3*sd(dataset_rmssdbase$resting_rmssd_log_cent))
#RMSSD Reactivity
dataset_rmssdbase$reactivity_rmssd_log[which(dataset_rmssdbase$reactivity_rmssd_log < (mean(dataset_rmssdbase$reactivity_rmssd_log) - 3*sd(dataset_rmssdbase$reactivity_rmssd_log)))] <- (mean(dataset_rmssdbase$reactivity_rmssd_log) - 3*sd(dataset_rmssdbase$reactivity_rmssd_log))
dataset_rmssdbase$reactivity_rmssd_log[which(dataset_rmssdbase$reactivity_rmssd_log > (mean(dataset_rmssdbase$reactivity_rmssd_log) + 3*sd(dataset_rmssdbase$reactivity_rmssd_log)))] <- (mean(dataset_rmssdbase$reactivity_rmssd_log) + 3*sd(dataset_rmssdbase$reactivity_rmssd_log))
#Reactivity MET
dataset_rmssdbase$reactivity_met_cent[which(dataset_rmssdbase$reactivity_met_cent < (mean(dataset_rmssdbase$reactivity_met_cent) - 3*sd(dataset_rmssdbase$reactivity_met_cent)))] <- (mean(dataset_rmssdbase$reactivity_met_cent) - 3*sd(dataset_rmssdbase$reactivity_met_cent))
dataset_rmssdbase$reactivity_met_cent[which(dataset_rmssdbase$reactivity_met_cent > (mean(dataset_rmssdbase$reactivity_met_cent) + 3*sd(dataset_rmssdbase$reactivity_met_cent)))] <- (mean(dataset_rmssdbase$reactivity_met_cent) + 3*sd(dataset_rmssdbase$reactivity_met_cent))
#Resting MET
dataset_rmssdbase$resting_met_cent[which(dataset_rmssdbase$resting_met_cent < (mean(dataset_rmssdbase$resting_met_cent) - 3*sd(dataset_rmssdbase$resting_met_cent)))] <- (mean(dataset_rmssdbase$resting_met_cent) - 3*sd(dataset_rmssdbase$resting_met_cent))
dataset_rmssdbase$resting_met_cent[which(dataset_rmssdbase$resting_met_cent > (mean(dataset_rmssdbase$resting_met_cent) + 3*sd(dataset_rmssdbase$resting_met_cent)))] <- (mean(dataset_rmssdbase$resting_met_cent) + 3*sd(dataset_rmssdbase$resting_met_cent))

#################################################################################################################
#base model RMSSD
########################################################################################################################

model_base <- lmerTest::lmer(rmssd_reactivity ~ resting_rmssd_log_cent + resting_met_cent + reactivity_met_cent + resting_rmssd_log_cent*reactivity_met_cent + conttime_cent + (1|id),data = dataset_rmssdbase)
tab_model(model_base, show.se = TRUE)

#################################################################################################################
#prepare dataset for covariates model
########################################################################################################################

dataset_rmssdco <- dataset[ , c("id", "stressmarker", "resting_rmssd", "reactivity_rmssd", "resting_met", "reactivity_met", "resting_supine", "reactivity_supine", "mvpa", "bmi", "chronicstress", "conttime", "caffeine_full", "nicotine_full", "alcohol_full", "ambientnoise_full", "anticipatedstress_full")]
dataset_rmssdco <- dataset_rmssdco[complete.cases(dataset_rmssdco),]

#RMSSD reactivity
dataset_rmssdco$rmssd_reactivity <- (dataset_rmssdco$reactivity_rmssd - dataset_rmssdco$resting_rmssd)

#Log transform RMSSD parameters
#RMSSD reactivity
dataset_rmssdco$reactivity_rmssd_log <- log(dataset_rmssdco$reactivity_rmssd)
#RMSSD resting
dataset_rmssdco$resting_rmssd_log <- log(dataset_rmssdco$resting_rmssd)

#supine position in minutes
dataset_rmssdco$resting_supine_minutes <- (dataset_rmssdco$resting_supine)*5
dataset_rmssdco$reactivity_supine_minutes <- (dataset_rmssdco$reactivity_supine)*5

#center independent variables
#Resting MET
dataset_rmssdco$resting_met_cent <- scale(dataset_rmssdco$resting_met, scale = FALSE)
#MVPA
dataset_rmssdco$mvpa_cent <- scale(dataset_rmssdco$mvpa, scale = FALSE)
#BMI
dataset_rmssdco$bmi_cent <- scale(dataset_rmssdco$bmi, scale = FALSE)
#Chronicstress
dataset_rmssdco$chronicstress_cent <- scale(dataset_rmssdco$chronicstress, scale = TRUE)
#RMSSD Resting
dataset_rmssdco$resting_rmssd_log_cent <- center(dataset_rmssdco$resting_rmssd_log, type = c("CWC"), cluster = dataset_rmssdco$id, value = NULL, as.na = NULL,
                                                   check = TRUE)
#Reactivity MET
dataset_rmssdco$reactivity_met_cent <- scale(dataset_rmssdco$reactivity_met, scale = FALSE)
#Conttime_hour
dataset_rmssdco$conttime_cent <- scale(dataset_rmssdco$conttime, scale = FALSE)

#outliers
#Reactivity RMSSD
dataset_rmssdco$rmssd_reactivity[which(dataset_rmssdco$rmssd_reactivity < (mean(dataset_rmssdco$rmssd_reactivity) - 3*sd(dataset_rmssdco$rmssd_reactivity)))] <- (mean(dataset_rmssdco$rmssd_reactivity) - 3*sd(dataset_rmssdco$rmssd_reactivity))
dataset_rmssdco$rmssd_reactivity[which(dataset_rmssdco$rmssd_reactivity > (mean(dataset_rmssdco$rmssd_reactivity) + 3*sd(dataset_rmssdco$rmssd_reactivity)))] <- (mean(dataset_rmssdco$rmssd_reactivity) + 3*sd(dataset_rmssdco$rmssd_reactivity))
#RMSSD resting
dataset_rmssdco$resting_rmssd_log_cent[which(dataset_rmssdco$resting_rmssd_log_cent < (mean(dataset_rmssdco$resting_rmssd_log_cent) - 3*sd(dataset_rmssdco$resting_rmssd_log_cent)))] <- (mean(dataset_rmssdco$resting_rmssd_log_cent) - 3*sd(dataset_rmssdco$resting_rmssd_log_cent))
dataset_rmssdco$resting_rmssd_log_cent[which(dataset_rmssdco$resting_rmssd_log_cent > (mean(dataset_rmssdco$resting_rmssd_log_cent) + 3*sd(dataset_rmssdco$resting_rmssd_log_cent)))] <- (mean(dataset_rmssdco$resting_rmssd_log_cent) + 3*sd(dataset_rmssdco$resting_rmssd_log_cent))
#RMSSD reactivity
dataset_rmssdco$reactivity_rmssd_log[which(dataset_rmssdco$reactivity_rmssd_log < (mean(dataset_rmssdco$reactivity_rmssd_log) - 3*sd(dataset_rmssdco$reactivity_rmssd_log)))] <- (mean(dataset_rmssdco$reactivity_rmssd_log) - 3*sd(dataset_rmssdco$reactivity_rmssd_log))
dataset_rmssdco$reactivity_rmssd_log[which(dataset_rmssdco$reactivity_rmssd_log > (mean(dataset_rmssdco$reactivity_rmssd_log) + 3*sd(dataset_rmssdco$reactivity_rmssd_log)))] <- (mean(dataset_rmssdco$reactivity_rmssd_log) + 3*sd(dataset_rmssdco$reactivity_rmssd_log))
#Reactivity MET
dataset_rmssdco$reactivity_met_cent[which(dataset_rmssdco$reactivity_met_cent < (mean(dataset_rmssdco$reactivity_met_cent) - 3*sd(dataset_rmssdco$reactivity_met_cent)))] <- (mean(dataset_rmssdco$reactivity_met_cent) - 3*sd(dataset_rmssdco$reactivity_met_cent))
dataset_rmssdco$reactivity_met_cent[which(dataset_rmssdco$reactivity_met_cent > (mean(dataset_rmssdco$reactivity_met_cent) + 3*sd(dataset_rmssdco$reactivity_met_cent)))] <- (mean(dataset_rmssdco$reactivity_met_cent) + 3*sd(dataset_rmssdco$reactivity_met_cent))
#Resting MET
dataset_rmssdco$resting_met_cent[which(dataset_rmssdco$resting_met_cent < (mean(dataset_rmssdco$resting_met_cent) - 3*sd(dataset_rmssdco$resting_met_cent)))] <- (mean(dataset_rmssdco$resting_met_cent) - 3*sd(dataset_rmssdco$resting_met_cent))
dataset_rmssdco$resting_met_cent[which(dataset_rmssdco$resting_met_cent > (mean(dataset_rmssdco$resting_met_cent) + 3*sd(dataset_rmssdco$resting_met_cent)))] <- (mean(dataset_rmssdco$resting_met_cent) + 3*sd(dataset_rmssdco$resting_met_cent))
#Resting Supine
dataset_rmssdco$resting_supine_minutes[which(dataset_rmssdco$resting_supine_minutes < (mean(dataset_rmssdco$resting_supine_minutes) - 3*sd(dataset_rmssdco$resting_supine_minutes)))] <- (mean(dataset_rmssdco$resting_supine_minutes) - 3*sd(dataset_rmssdco$resting_supine_minutes))
dataset_rmssdco$resting_supine_minutes[which(dataset_rmssdco$resting_supine_minutes > (mean(dataset_rmssdco$resting_supine_minutes) + 3*sd(dataset_rmssdco$resting_supine_minutes)))] <- (mean(dataset_rmssdco$resting_supine_minutes) + 3*sd(dataset_rmssdco$resting_supine_minutes))
#Reactivity Supine
dataset_rmssdco$reactivity_supine_minutes[which(dataset_rmssdco$reactivity_supine_minutes < (mean(dataset_rmssdco$reactivity_supine_minutes) - 3*sd(dataset_rmssdco$reactivity_supine_minutes)))] <- (mean(dataset_rmssdco$reactivity_supine_minutes) - 3*sd(dataset_rmssdco$reactivity_supine_minutes))
dataset_rmssdco$reactivity_supine_minutes[which(dataset_rmssdco$reactivity_supine_minutes > (mean(dataset_rmssdco$reactivity_supine_minutes) + 3*sd(dataset_rmssdco$reactivity_supine_minutes)))] <- (mean(dataset_rmssdco$reactivity_supine_minutes) + 3*sd(dataset_rmssdco$reactivity_supine_minutes))
#MVPA
dataset_rmssdco$mvpa_cent[which(dataset_rmssdco$mvpa_cent < (mean(dataset_rmssdco$mvpa_cent) - 3*sd(dataset_rmssdco$mvpa_cent)))] <- (mean(dataset_rmssdco$mvpa_cent) - 3*sd(dataset_rmssdco$mvpa_cent))
dataset_rmssdco$mvpa_cent[which(dataset_rmssdco$mvpa_cent > (mean(dataset_rmssdco$mvpa_cent) + 3*sd(dataset_rmssdco$mvpa_cent)))] <- (mean(dataset_rmssdco$mvpa_cent) + 3*sd(dataset_rmssdco$mvpa_cent))
#BMI
dataset_rmssdco$bmi_cent[which(dataset_rmssdco$bmi_cent < (mean(dataset_rmssdco$bmi_cent) - 3*sd(dataset_rmssdco$bmi_cent)))] <- (mean(dataset_rmssdco$bmi_cent) - 3*sd(dataset_rmssdco$bmi_cent))
dataset_rmssdco$bmi_cent[which(dataset_rmssdco$bmi_cent > (mean(dataset_rmssdco$bmi_cent) + 3*sd(dataset_rmssdco$bmi_cent)))] <- (mean(dataset_rmssdco$bmi_cent) + 3*sd(dataset_rmssdco$bmi_cent))
#Chronicstress
dataset_rmssdco$chronicstress_cent[which(dataset_rmssdco$chronicstress_cent < (mean(dataset_rmssdco$chronicstress_cent) - 3*sd(dataset_rmssdco$chronicstress_cent)))] <- (mean(dataset_rmssdco$chronicstress_cent) - 3*sd(dataset_rmssdco$chronicstress_cent))
dataset_rmssdco$chronicstress_cent[which(dataset_rmssdco$chronicstress_cent > (mean(dataset_rmssdco$chronicstress_cent) + 3*sd(dataset_rmssdco$chronicstress_cent)))] <- (mean(dataset_rmssdco$chronicstress_cent) + 3*sd(dataset_rmssdco$chronicstress_cent))

#rescale
dataset_rmssdco$conttime_cent <- dataset_rmssdco$conttime_cent/24
dataset_rmssdco$mvpa_cent <- dataset_rmssdco$mvpa_cent/60

################################################################################################################
#Models with covariates
########################################################################################################################

model_covariates <- lmerTest::lmer(rmssd_reactivity ~ resting_rmssd_log_cent + resting_met_cent + reactivity_met_cent + resting_rmssd_log_cent*reactivity_met_cent + resting_supine_minutes + reactivity_supine_minutes + conttime_cent + caffeine_full + alcohol_full + anticipatedstress_full + ambientnoise_full + mvpa_cent + bmi_cent + chronicstress_cent + (1|id),data = dataset_rmssdco)
tab_model(model_covariates, show.se = TRUE)

################################################################################################################
#Figures
########################################################################################################################

dataset_rmssdco$predicted <- predict(model_covariates)

tiff("Figure2__new.tiff", units="in", width=5, height=5, res=300)
ggplot(data=dataset_rmssdco, aes(x=resting_rmssd_log_cent, y=predicted, group=factor(id), color="gray"), legend=FALSE) +
  geom_smooth(method=lm, se=FALSE, fullrange=TRUE, lty=1, size=.5, color="grey40") +
  geom_smooth(aes(group=1), method=lm, se=FALSE, fullrange=FALSE, lty=1, size=2, color="blue") +
  xlab("Resting lnRMSSD") + ylab("Predicted RMSSD Reactivity (ms)") +
  theme_classic() +
  theme(axis.title=element_text(size=18),
        axis.text=element_text(size=18),
        plot.title=element_text(size=18, hjust=.5)) +
  ggtitle("Within-Person Association Plot\nResting HRV & HRV Reactivity")
dev.off()
